03 analysis Imbizo 2023
Introduction¶
General¶
- Please read this introductory notebook to get familiar with the Brainwide dataset and the task structure.
Once this is done, please read the specific instructions below.
What are the data of interest¶
Choice 💭¶
In this notebook, we will teach you how to run analysis, that will revolve around a variable of interest: the choice.
The choice is the decision the mouse makes to turn the wheel clock-wise (CW) or counter-clock-wise (CCW), in order to correctly bring a visual stimulus presented on the left or on the right to the center of the screen respectively.
The choice values are encoded and represented as:
- value
+1= CW (right), painted in blue - value
-1= CCW (left), painted in red
Wheel ⭕¶
As we will want to predict the choice of the animal, we need to chose a time window that is not contaminated by movement. To check for such lack of movement, you will explore the wheel data, and plot its speed and position.
First movement 🔽¶
The time window for analysis will be taken before the turning of the wheel has started (as such turning is in itself indicative of the choice) . You will use the onset of the first movement as detected on the wheel trace to define your time window of interest.
Spikes and clusters 🧠¶
You will use the spikes from selected single-units (clusters). We will call these units throughout the text.
What questions will each analysis answer¶
All of the analysis below use the neural activity in a time window of ~100ms before the first movement.
Encoding 🔛¶
If you were to rely on a single unit to predict the choice value, which one would you chose, and what would be the magnitude of the difference between its response to left or right choices?
- Note: This analysis concatenates the activity across trials
Manifold 📈¶
Using all units in a given brain region, when exactly would you be able to differentiate between the left or right choice in our time window of interest?
- Note: This analysis concatenates the activity across trials
Decoding ⏭¶
For a given trial, if you were a downstream neuron receiving the activity from all units in a given region as input, would you be able to differentiate left versus right choices?
Explore further 🌍¶
We will give at the end a set of ideas for you to explore the data further.
Explore the data before you start¶
Important reference:
- You can see the data for the example recording through this visualisation website.
Browse through a few trial in the single trial overview panel, and notice the timing between the different events, notably between the stim on, first move and feedback.
Notice how the wheel position changes during the trial: it is still before stim on and moves between the first move and feedback.
Notice that the wheel trace is very different between left and right choices, and that the wheel is nearly not moving before the first movement onset (dashed line):
On the website, browse the response of single good unit - here shown for unit #55. You can see that some units are clearly responding selectively to a particular choice after the first movement onset.
Getting help¶
Please post questions on :
- Dandi / NWB help desk: https://github.com/dandi/helpdesk/discussions/
- IBL help desk: https://neurostars.org/tag/ibl
Installation¶
References:
- The IBL dataset can be found on the Dandi archive here.
- IBL documentation website
- BW main code repository: https://github.com/int-brain-lab/paper-brain-wide-map
In [7]:
Copied!
%%capture
!pip install tabulate
!pip install ibllib
!git clone -b imbizo https://github.com/int-brain-lab/paper-brain-wide-map.git
%cd ./paper-brain-wide-map
!pip install -e .
%%capture
!pip install tabulate
!pip install ibllib
!git clone -b imbizo https://github.com/int-brain-lab/paper-brain-wide-map.git
%cd ./paper-brain-wide-map
!pip install -e .
Code init and data selection¶
Imports and inits¶
In [3]:
Copied!
# ---------------------------------------------------
# Import
# DANDI
import fsspec
import h5py
from fsspec.implementations.cached import CachingFileSystem
from pynwb import NWBHDF5IO
from dandi.dandiapi import DandiAPIClient
# IBL
from brainbox.population.decode import get_spike_counts_in_bins
from brainbox.ephys_plots import plot_brain_regions
from brainbox.behavior.wheel import interpolate_position, velocity_filtered
from brainbox.task.trials import get_event_aligned_raster, get_psth
from ibllib.atlas import AllenAtlas
from brainwidemap import bwm_query
from brainwidemap.imbizo.encoding_functions import get_choice_time_shuffle
# Generic
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
# Init
ba = AllenAtlas()
# ---------------------------------------------------
# Import
# DANDI
import fsspec
import h5py
from fsspec.implementations.cached import CachingFileSystem
from pynwb import NWBHDF5IO
from dandi.dandiapi import DandiAPIClient
# IBL
from brainbox.population.decode import get_spike_counts_in_bins
from brainbox.ephys_plots import plot_brain_regions
from brainbox.behavior.wheel import interpolate_position, velocity_filtered
from brainbox.task.trials import get_event_aligned_raster, get_psth
from ibllib.atlas import AllenAtlas
from brainwidemap import bwm_query
from brainwidemap.imbizo.encoding_functions import get_choice_time_shuffle
# Generic
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
# Init
ba = AllenAtlas()
Downloading: /home/jovyan/Downloads/ONE/openalyx.internationalbrainlab.org/histology/ATLAS/Needles/Allen/average_template_25.nrrd Bytes: 32998960
100%|██████████| 31.470260620117188/31.470260620117188 [00:04<00:00, 7.51it/s]
Downloading: /home/jovyan/Downloads/ONE/openalyx.internationalbrainlab.org/histology/ATLAS/Needles/Allen/annotation_25.nrrd Bytes: 4035363
100%|██████████| 3.848422050476074/3.848422050476074 [00:00<00:00, 4.45it/s]
2023-09-04 21:57:51.430 INFO [atlas.py:1335] Computing brain atlas annotations lookup table 2023-09-04 21:58:02.152 INFO [atlas.py:1359] Cached remapping file /home/jovyan/Downloads/ONE/openalyx.internationalbrainlab.org/histology/ATLAS/Needles/Allen/annotation_25_lut_v01.npz ...
Load the NWB file for a single session¶
🔀 TODO How to search on DANDI for available sessions.
Define a session filepath of interest, and read the associated NWB file :
In [5]:
Copied!
# get S3 URL for NWB file on DANDI
dandiset_id = '000409'
filepath = 'sub-KS042/sub-KS042_ses-07dc4b76-5b93-4a03-82a0-b3d9cc73f412_behavior+ecephys+image.nwb'
cache_dir = "/tmp/fsspec_cache" # Local folder for the cache
with DandiAPIClient() as client:
asset = client.get_dandiset(dandiset_id, 'draft').get_asset_by_path(filepath)
s3_url = asset.get_content_url(follow_redirects=1, strip_query=True)
# configure CachingFileSystem for streaming
cfs = CachingFileSystem(
fs=fsspec.filesystem("http"),
cache_storage=cache_dir,
)
file_system = cfs.open(s3_url, "rb")
file = h5py.File(file_system)
io = NWBHDF5IO(file=file, load_namespaces=True)
# Read in the NWB file
nwb = io.read()
# get S3 URL for NWB file on DANDI
dandiset_id = '000409'
filepath = 'sub-KS042/sub-KS042_ses-07dc4b76-5b93-4a03-82a0-b3d9cc73f412_behavior+ecephys+image.nwb'
cache_dir = "/tmp/fsspec_cache" # Local folder for the cache
with DandiAPIClient() as client:
asset = client.get_dandiset(dandiset_id, 'draft').get_asset_by_path(filepath)
s3_url = asset.get_content_url(follow_redirects=1, strip_query=True)
# configure CachingFileSystem for streaming
cfs = CachingFileSystem(
fs=fsspec.filesystem("http"),
cache_storage=cache_dir,
)
file_system = cfs.open(s3_url, "rb")
file = h5py.File(file_system)
io = NWBHDF5IO(file=file, load_namespaces=True)
# Read in the NWB file
nwb = io.read()
/opt/conda/lib/python3.10/site-packages/hdmf/spec/namespace.py:531: UserWarning: Ignoring cached namespace 'core' version 2.5.0 because version 2.6.0-alpha is already loaded.
warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."
Format the data necessary for the analysis¶
Spikes¶
Spikes and clusters information are saved in a single table, where each row corresponds to a cluster, and spike times are contained within the column "spike_times".
In [8]:
Copied!
# ---------------------------------------------------
# Load spike data
units = nwb.units.to_dataframe()
# Show first 3 rows of dataframe
print(units.head(3).to_markdown())
# ---------------------------------------------------
# Load spike data
units = nwb.units.to_dataframe()
# Show first 3 rows of dataframe
print(units.head(3).to_markdown())
| id | spike_times | unit_name | presence_ratio_standard_deviation | contamination | noise_cutoff | mean_relative_depth | sliding_refractory_period_violation | cosmos_location | maximum_amplitude | maximum_amplitude_channel | alternative_contamination | probe_name | label | spike_relative_depths | beryl_location | median_amplitude | drift | minimum_amplitude | spike_count | firing_rate | missed_spikes_estimate | standard_deviation_amplitude | allen_location | spike_amplitudes | presence_ratio | |-----:|:-----------------------------------------------------------------|------------:|------------------------------------:|----------------:|---------------:|----------------------:|--------------------------------------:|:------------------|--------------------:|----------------------------:|----------------------------:|:-------------|---------:|:------------------------------------------------------------------------|:-----------------|-------------------:|----------:|--------------------:|--------------:|--------------:|-------------------------:|-------------------------------:|:-----------------|:------------------------------------------------------------------------|-----------------:| | 0 | [1.38666667e-02 4.66000000e-02 7.73000000e-02 ... 5.36146927e+03 | 52 | 0.974625 | 0 | 36.5833 | 40 | 0 | MB | 0.000164929 | 3 | 0 | probe01 | 0.333333 | [6.39112941e-05 8.90963566e-05 9.62116423e-05 ... 7.58570206e-05 | MRN | 0.000103959 | 2755.33 | 9.66621e-05 | 326 | 0.0608029 | nan | 0.92908 | MRN | [6.39112941e-05 8.90963566e-05 9.62116423e-05 ... 7.58570206e-05 | 0.392924 | | | 5.36150370e+03 5.36151427e+03] | | | | | | | | | | | | | 6.26894223e-05 7.55858358e-05] | | | | | | | | | | 6.26894223e-05 7.55858358e-05] | | | 1 | [1.48000000e-02 3.66000000e-02 6.96000000e-02 ... 5.32395790e+03 | 173 | 36.8049 | 0.102952 | 370.745 | 40 | 1 | MB | 0.000261255 | 2 | 0.0878241 | probe01 | 0.666667 | [5.89352876e-05 1.17742655e-04 1.03836896e-04 ... 5.56572846e-05 | MRN | 0.000115008 | 314121 | 0.00010324 | 27276 | 5.0873 | 0.5 | 1.03926 | MRN | [5.89352876e-05 1.17742655e-04 1.03836896e-04 ... 5.56572846e-05 | 0.968343 | | | 5.34243170e+03 5.34347823e+03] | | | | | | | | | | | | | 6.54943406e-05 5.56091879e-05] | | | | | | | | | | 6.54943406e-05 5.56091879e-05] | | | 2 | [1.48333333e-02 3.66333333e-02 6.96333333e-02 ... 5.36154460e+03 | 175 | 66.5869 | 0.0545487 | 74.5426 | 40 | 1 | MB | 0.0003058 | 2 | 0.0464068 | probe01 | 0.666667 | [0.00017621 0.00029232 0.00027048 ... 0.00030599 0.00027077 0.00027227] | MRN | 0.00013126 | 375654 | 9.6097e-05 | 32988 | 6.15266 | 0.239371 | 1.5242 | MRN | [0.00017621 0.00029232 0.00027048 ... 0.00030599 0.00027077 0.00027227] | 0.765363 | | | 5.36156417e+03 5.36158920e+03] | | | | | | | | | | | | | | | | | | | | | | | | |
/opt/conda/lib/python3.10/site-packages/tabulate/__init__.py:107: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison (len(row) >= 1 and row[0] == SEPARATING_LINE) /opt/conda/lib/python3.10/site-packages/tabulate/__init__.py:108: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison or (len(row) >= 2 and row[1] == SEPARATING_LINE)
Note that the clusters/spikes acquired on different Neuropixels probes in a session are concatenated in this table.
To differentiate between the probes, use the column probe_name.
In this example, there are 2 probes labelled probe00 and probe01:
In [9]:
Copied!
units.probe_name.unique()
units.probe_name.unique()
Out[9]:
array(['probe01', 'probe00'], dtype=object)
For this analysis, we will use the clusters from probe00, and restrict only to good clusters.
In [10]:
Copied!
# ---------------------------------------------------
# Restrict to only good clusters and one probe for this analysis
probe_label = 'probe00'
good_units = units[(units.label == 1) & (units.probe_name == probe_label)]
# We sort by unit number
good_units = good_units.sort_values('unit_name')
# Reset the index so we start from 0
good_units = good_units.reset_index()
# ---------------------------------------------------
# N neuronal units in total
num_neuron = len(good_units)
# Show example data for first 3 units
print(f'We will be working with {num_neuron} good units \n')
print(good_units.head(3).to_markdown())
# ---------------------------------------------------
# Restrict to only good clusters and one probe for this analysis
probe_label = 'probe00'
good_units = units[(units.label == 1) & (units.probe_name == probe_label)]
# We sort by unit number
good_units = good_units.sort_values('unit_name')
# Reset the index so we start from 0
good_units = good_units.reset_index()
# ---------------------------------------------------
# N neuronal units in total
num_neuron = len(good_units)
# Show example data for first 3 units
print(f'We will be working with {num_neuron} good units \n')
print(good_units.head(3).to_markdown())
We will be working with 51 good units | | id | spike_times | unit_name | presence_ratio_standard_deviation | contamination | noise_cutoff | mean_relative_depth | sliding_refractory_period_violation | cosmos_location | maximum_amplitude | maximum_amplitude_channel | alternative_contamination | probe_name | label | spike_relative_depths | beryl_location | median_amplitude | drift | minimum_amplitude | spike_count | firing_rate | missed_spikes_estimate | standard_deviation_amplitude | allen_location | spike_amplitudes | presence_ratio | |---:|-----:|:-----------------------------------------------------------------|------------:|------------------------------------:|----------------:|---------------:|----------------------:|--------------------------------------:|:------------------|--------------------:|----------------------------:|----------------------------:|:-------------|--------:|:------------------------------------------------------------------------|:-----------------|-------------------:|---------:|--------------------:|--------------:|--------------:|-------------------------:|-------------------------------:|:-----------------|:------------------------------------------------------------------------|-----------------:| | 0 | 695 | [2.21643333e+00 4.57046667e+00 1.55576000e+01 ... 5.33761487e+03 | 545 | 35.5377 | 1.17832 | -0.522168 | 2100 | 1 | MB | 0.000999434 | 209 | 0.661216 | probe00 | 1 | [0.00013078 0.00013274 0.00013505 ... 0.00013053 0.00013372 0.00015646] | SCm | 0.000692823 | 27423.7 | 0.00032281 | 23917 | 4.46082 | 0.000288285 | 1.6835 | SCiw | [0.00013078 0.00013274 0.00013505 ... 0.00013053 0.00013372 0.00015646] | 0.988827 | | | | 5.35226167e+03 5.35245793e+03] | | | | | | | | | | | | | | | | | | | | | | | | | | 1 | 689 | [1.55826667e+00 3.84840000e+00 6.48830000e+00 1.10399000e+01 | 548 | 96.6628 | 0.0478893 | -0.736521 | 2060 | 1 | MB | 0.000419182 | 205 | 0.042819 | probe00 | 1 | [0.00017287 0.00016987 0.00016825 0.00018331 0.00017805 0.00016998 | SCm | 0.000244842 | 593508 | 0.00011452 | 30978 | 5.77778 | 0.00139621 | 1.3312 | SCiw | [0.00017287 0.00016987 0.00016825 0.00018331 0.00017805 0.00016998 | 0.873371 | | | | 1.12341333e+01 1.30420333e+01 1.30847667e+01 1.35472667e+01 | | | | | | | | | | | | | 0.00017017 0.00018654 0.00017474 0.00016566 0.00019011 0.00016557 | | | | | | | | | | 0.00017017 0.00018654 0.00017474 0.00016566 0.00019011 0.00016557 | | | | | 1.37510000e+01 1.51390000e+01 1.72040333e+01 1.80920333e+01 | | | | | | | | | | | | | 0.00020082 0.00021913 0.00018343 0.00018929 0.00017776 0.00019305 | | | | | | | | | | 0.00020082 0.00021913 0.00018343 0.00018929 0.00017776 0.00019305 | | | | | 1.85952667e+01 1.94287000e+01 2.10797000e+01 2.15468333e+01 | | | | | | | | | | | | | 0.00021832 0.0001679 0.00016918 0.00019594 0.00017234 0.00017121 | | | | | | | | | | 0.00021832 0.0001679 0.00016918 0.00019594 0.00017234 0.00017121 | | | | | 2.41510667e+01 2.58954000e+01 2.82780000e+01 2.92861333e+01 | | | | | | | | | | | | | 0.00017205 0.00018681 0.00016775 0.00018708 0.00017999 0.00021725 | | | | | | | | | | 0.00017205 0.00018681 0.00016775 0.00018708 0.00017999 0.00021725 | | | | | 3.03733333e+01 3.08643333e+01 3.16312333e+01 3.16655667e+01 | | | | | | | | | | | | | 0.00020638 0.00016855 0.00017333 0.00018273 0.00017783 0.00018538 | | | | | | | | | | 0.00020638 0.00016855 0.00017333 0.00018273 0.00017783 0.00018538 | | | | | 3.18054333e+01 3.38607000e+01 3.50156333e+01 3.79925000e+01 | | | | | | | | | | | | | 0.00019111 0.00020817 0.00017819 0.00020045 0.00020139 0.00016832 | | | | | | | | | | 0.00019111 0.00020817 0.00017819 0.00020045 0.00020139 0.00016832 | | | | | 3.87990000e+01 3.88412333e+01 4.10038000e+01 4.16848667e+01 | | | | | | | | | | | | | 0.00017348 0.00020975 0.00020225 0.00020365 0.00019568 0.00017773 | | | | | | | | | | 0.00017348 0.00020975 0.00020225 0.00020365 0.00019568 0.00017773 | | | | | 4.41783333e+01 4.41813667e+01 4.43800667e+01 4.67868000e+01 | | | | | | | | | | | | | 0.00016582 0.00018897 0.00021721 0.00017155 0.00020643 0.00017612 | | | | | | | | | | 0.00016582 0.00018897 0.00021721 0.00017155 0.00020643 0.00017612 | | | | | 4.74190667e+01 4.86860667e+01 4.92406000e+01 5.02259667e+01 | | | | | | | | | | | | | 0.00019797 0.00016867 0.00017643 0.00017057 0.00017489 0.00020119 | | | | | | | | | | 0.00019797 0.00016867 0.00017643 0.00017057 0.00017489 0.00020119 | | | | | 5.10363667e+01 5.40579000e+01 5.78586667e+01 5.81384000e+01 | | | | | | | | | | | | | 0.00018548 0.00019693 0.0002232 0.00017938 0.00023602 0.00016944 | | | | | | | | | | 0.00018548 0.00019693 0.0002232 0.00017938 0.00023602 0.00016944 | | | | | 6.01455000e+01 6.08485000e+01 6.42386333e+01 6.95582000e+01 | | | | | | | | | | | | | 0.00022088 0.00016882 0.00018777 0.00017846 0.00018092 0.00019765 | | | | | | | | | | 0.00022088 0.00016882 0.00018777 0.00017846 0.00018092 0.00019765 | | | | | 7.40021667e+01 7.82306000e+01 8.12250333e+01 8.22878000e+01 | | | | | | | | | | | | | 0.00022248 0.00018469 0.0002265 0.00017194 0.00017227 0.00017035 | | | | | | | | | | 0.00022248 0.00018469 0.0002265 0.00017194 0.00017227 0.00017035 | | | | | 8.66895000e+01 9.78342333e+01 9.86141000e+01 9.91412667e+01 | | | | | | | | | | | | | 0.00017837 0.00021148 0.00022106 0.00020092 0.00020107 0.00021783 | | | | | | | | | | 0.00017837 0.00021148 0.00022106 0.00020092 0.00020107 0.00021783 | | | | | 1.00035767e+02 1.00653800e+02 1.01900267e+02 1.03160367e+02 | | | | | | | | | | | | | 0.00020034 0.00017961 0.00017473 0.00022378 0.00017427 0.00019447 | | | | | | | | | | 0.00020034 0.00017961 0.00017473 0.00022378 0.00017427 0.00019447 | | | | | 1.06653300e+02 1.09354567e+02 1.12514300e+02 1.12569400e+02 | | | | | | | | | | | | | 0.00018281 0.00026028 0.00017552 0.00021767 0.00018056 0.00017286 | | | | | | | | | | 0.00018281 0.00026028 0.00017552 0.00021767 0.00018056 0.00017286 | | | | | 1.12622333e+02 1.12714800e+02 1.12856167e+02 1.13374700e+02 | | | | | | | | | | | | | 0.00021191 0.00020853 0.00020159 0.0001903 0.00019292 0.00016989 | | | | | | | | | | 0.00021191 0.00020853 0.00020159 0.0001903 0.00019292 0.00016989 | | | | | 1.13698367e+02 1.15029400e+02 1.19231267e+02 1.21053133e+02 | | | | | | | | | | | | | 0.00017524 0.00020974 0.00018376 0.00020333 0.00021557 0.00018139 | | | | | | | | | | 0.00017524 0.00020974 0.00018376 0.00020333 0.00021557 0.00018139 | | | | | 1.21972000e+02 1.24376633e+02 1.24484033e+02 1.24989867e+02 | | | | | | | | | | | | | 0.00020685 0.00021795 0.00028413 0.00016976 0.00021461 0.00023207 | | | | | | | | | | 0.00020685 0.00021795 0.00028413 0.00016976 0.00021461 0.00023207 | | | | | 1.27305467e+02 1.30272800e+02 1.30655533e+02 1.30739500e+02 | | | | | | | | | | | | | 0.00017286 0.000201 0.0002164 0.00018295 0.00022594 0.00018511 | | | | | | | | | | 0.00017286 0.000201 0.0002164 0.00018295 0.00022594 0.00018511 | | | | | 1.30791000e+02 1.30939733e+02 1.33360900e+02 1.33370100e+02 | | | | | | | | | | | | | 0.00020801 0.00016879 0.00017408 0.00018744 0.00022867 0.00018829 | | | | | | | | | | 0.00020801 0.00016879 0.00017408 0.00018744 0.00022867 0.00018829 | | | | | 1.36410567e+02 1.36720967e+02 1.37729100e+02 1.38980800e+02 | | | | | | | | | | | | | 0.00023807 0.00020105 0.0001906 0.00017126 0.00016927 0.00020795 | | | | | | | | | | 0.00023807 0.00020105 0.0001906 0.00017126 0.00016927 0.00020795 | | | | | 1.38998200e+02 1.39344733e+02 1.44716367e+02 1.44726900e+02 | | | | | | | | | | | | | 0.00021704 0.0002175 0.00030227 0.00017033 0.00021817 0.00021402 | | | | | | | | | | 0.00021704 0.0002175 0.00030227 0.00017033 0.00021817 0.00021402 | | | | | 1.47877400e+02 1.48346667e+02 1.49946667e+02 1.49981033e+02 | | | | | | | | | | | | | 0.0002373 0.00018441 0.00023337 0.0001814 0.00020103 0.00019011 | | | | | | | | | | 0.0002373 0.00018441 0.00023337 0.0001814 0.00020103 0.00019011 | | | | | 1.51224900e+02 1.51905267e+02 1.52025267e+02 1.53933300e+02 | | | | | | | | | | | | | 0.00022198 0.00020894 0.00025814 0.00017253 0.00016738 0.00020995 | | | | | | | | | | 0.00022198 0.00020894 0.00025814 0.00017253 0.00016738 0.00020995 | | | | | 1.55030100e+02 1.59071833e+02 1.60224700e+02 1.60430900e+02 | | | | | | | | | | | | | 0.00020144 0.00021374 0.0002009 0.00016919 0.0002001 0.00018027 | | | | | | | | | | 0.00020144 0.00021374 0.0002009 0.00016919 0.0002001 0.00018027 | | | | | 1.67128000e+02 1.67969200e+02 1.69069500e+02 1.69152600e+02 | | | | | | | | | | | | | 0.00020637 0.00018815 0.00017311 0.00017186 0.00017976 0.00018765 | | | | | | | | | | 0.00020637 0.00018815 0.00017311 0.00017186 0.00017976 0.00018765 | | | | | 1.69264933e+02 1.71874500e+02 1.73502933e+02 1.76165200e+02 | | | | | | | | | | | | | 0.00018957 0.00018183 0.00017148 0.00024981 0.00017877 0.00016848 | | | | | | | | | | 0.00018957 0.00018183 0.00017148 0.00024981 0.00017877 0.00016848 | | | | | 1.76485500e+02 1.79894667e+02 1.80402067e+02 1.88235333e+02 | | | | | | | | | | | | | 0.00021596 0.00024345 0.00017291 0.00019193 0.00030302 0.00022514 | | | | | | | | | | 0.00021596 0.00024345 0.00017291 0.00019193 0.00030302 0.00022514 | | | | | 1.97494900e+02 1.99366800e+02 1.99601033e+02 2.01365333e+02 | | | | | | | | | | | | | 0.00024498 0.00022841 0.00016673 0.00016782 0.00022167 0.00019854 | | | | | | | | | | 0.00024498 0.00022841 0.00016673 0.00016782 0.00022167 0.00019854 | | | | | 2.03367367e+02 2.33280667e+02 2.35431800e+02 2.52358633e+02 | | | | | | | | | | | | | 0.00026057 0.00016813 0.00019533 0.0002019 0.0001669 0.00020812 | | | | | | | | | | 0.00026057 0.00016813 0.00019533 0.0002019 0.0001669 0.00020812 | | | | | 2.55002633e+02 2.55079733e+02 2.57506867e+02 2.57583267e+02 | | | | | | | | | | | | | 0.00021787 0.00023673 0.00017533 0.0002322 0.00024007 0.00019833 | | | | | | | | | | 0.00021787 0.00023673 0.00017533 0.0002322 0.00024007 0.00019833 | | | | | 2.65569467e+02 2.65798667e+02 2.69909533e+02 3.08805900e+02 | | | | | | | | | | | | | 0.00016675 0.00019151 0.00017359 0.00023753 0.00018435 0.00019643 | | | | | | | | | | 0.00016675 0.00019151 0.00017359 0.00023753 0.00018435 0.00019643 | | | | | 3.20413067e+02 3.32747967e+02 3.40795800e+02 3.51130000e+02 | | | | | | | | | | | | | 0.00016636 0.0001859 0.0001771 0.00019928 0.00018085 0.00020749 | | | | | | | | | | 0.00016636 0.0001859 0.0001771 0.00019928 0.00018085 0.00020749 | | | | | 3.59876333e+02 3.62159800e+02 3.62355167e+02 3.72977767e+02 | | | | | | | | | | | | | 0.00018921 0.00020761 0.0002393 0.00026477 0.00016567 0.00016974 | | | | | | | | | | 0.00018921 0.00020761 0.0002393 0.00026477 0.00016567 0.00016974 | | | | | 4.33538533e+02 4.33727033e+02 4.63451267e+02 4.94539467e+02 | | | | | | | | | | | | | 0.00019104 0.00016896 0.00020921 0.00021288 0.00016874 0.00020456 | | | | | | | | | | 0.00019104 0.00016896 0.00020921 0.00021288 0.00016874 0.00020456 | | | | | 5.22367667e+02 5.24235733e+02 5.26998033e+02 5.30443967e+02 | | | | | | | | | | | | | 0.00017128 0.00018062 0.00016673 0.00020854 0.0002254 0.00018631 | | | | | | | | | | 0.00017128 0.00018062 0.00016673 0.00020854 0.0002254 0.00018631 | | | | | 5.83249767e+02 6.51383567e+02 7.37111800e+02 1.20915713e+03 | | | | | | | | | | | | | 0.00016656 0.00021065 0.00018781 0.00022755 0.00019399 0.0002304 | | | | | | | | | | 0.00016656 0.00021065 0.00018781 0.00022755 0.00019399 0.0002304 | | | | | 1.23220240e+03 1.32069387e+03 1.32469273e+03 1.47550400e+03 | | | | | | | | | | | | | 0.00022553 0.00020407 0.0001851 0.00017808 0.00019342 0.00017096 | | | | | | | | | | 0.00022553 0.00020407 0.0001851 0.00017808 0.00019342 0.00017096 | | | | | 1.54771147e+03 1.57636140e+03 1.57647637e+03 1.59130867e+03 | | | | | | | | | | | | | 0.00020796 0.00016656 0.00019078 0.00022391 0.00021116 0.00017101 | | | | | | | | | | 0.00020796 0.00016656 0.00019078 0.00022391 0.00021116 0.00017101 | | | | | 1.59400243e+03 1.64592590e+03 1.76955750e+03 1.76972630e+03 | | | | | | | | | | | | | 0.00017826 0.00021425 0.00019773 0.00017841 0.00017222 0.00021095 | | | | | | | | | | 0.00017826 0.00021425 0.00019773 0.00017841 0.00017222 0.00021095 | | | | | 1.87403363e+03 1.93367053e+03 1.96506923e+03 2.14627490e+03 | | | | | | | | | | | | | 0.00017957 0.00017634 0.00016557 0.00019645 0.00017277 0.00016716 | | | | | | | | | | 0.00017957 0.00017634 0.00016557 0.00019645 0.00017277 0.00016716 | | | | | 2.15130203e+03 2.17374433e+03 2.18568800e+03 2.19219960e+03 | | | | | | | | | | | | | 0.0001753 0.00017198 0.00022106 0.00019481 0.00020494 0.00018171 | | | | | | | | | | 0.0001753 0.00017198 0.00022106 0.00019481 0.00020494 0.00018171 | | | | | 2.21064323e+03 2.21332187e+03 2.25844063e+03 2.39620050e+03 | | | | | | | | | | | | | 0.00019048 0.00017358 0.0001664 0.00017013 0.00021874 0.00020574 | | | | | | | | | | 0.00019048 0.00017358 0.0001664 0.00017013 0.00021874 0.00020574 | | | | | 2.47417863e+03 2.47459017e+03 2.48185903e+03 2.55221283e+03 | | | | | | | | | | | | | 0.00017252 0.00018783 0.00017031 0.0001775 0.00018547 0.00018115 | | | | | | | | | | 0.00017252 0.00018783 0.00017031 0.0001775 0.00018547 0.00018115 | | | | | 2.61064567e+03 2.66344120e+03 2.73972360e+03 2.78752057e+03 | | | | | | | | | | | | | 0.00017571 0.00017215 0.00021133 0.00017497 0.00017999 0.00019389 | | | | | | | | | | 0.00017571 0.00017215 0.00021133 0.00017497 0.00017999 0.00019389 | | | | | 2.79321597e+03 2.81803827e+03 2.87586497e+03 2.96569150e+03 | | | | | | | | | | | | | 0.00017093 0.000174 0.000208 0.00020423 0.0001751 0.00018906 | | | | | | | | | | 0.00017093 0.000174 0.000208 0.00020423 0.0001751 0.00018906 | | | | | 2.98025360e+03 2.98679530e+03 2.98998490e+03 2.99092087e+03 | | | | | | | | | | | | | 0.00017228 0.00018614 0.00019254 0.00019652 0.00016597 0.00016919 | | | | | | | | | | 0.00017228 0.00018614 0.00019254 0.00019652 0.00016597 0.00016919 | | | | | 2.99547047e+03 3.01261177e+03 3.04886130e+03 3.05935790e+03 | | | | | | | | | | | | | 0.0001866 0.00025445 0.00018718 0.00020683 0.00020531 0.00018404 | | | | | | | | | | 0.0001866 0.00025445 0.00018718 0.00020683 0.00020531 0.00018404 | | | | | 3.05935957e+03 3.06440887e+03 3.06451397e+03 3.07541983e+03 | | | | | | | | | | | | | 0.00020755 0.00016832 0.00018836 0.00017131 0.00019283 0.000169 | | | | | | | | | | 0.00020755 0.00016832 0.00018836 0.00017131 0.00019283 0.000169 | | | | | 3.07862420e+03 3.08183917e+03 3.09398460e+03 3.10382020e+03 | | | | | | | | | | | | | 0.00016738 0.00018245 0.00018232 0.00018056 0.00021608 0.00020182 | | | | | | | | | | 0.00016738 0.00018245 0.00018232 0.00018056 0.00021608 0.00020182 | | | | | 3.13272837e+03 3.17651267e+03 3.23143157e+03 3.24501430e+03 | | | | | | | | | | | | | 0.00020881 0.00017304 0.00018072 0.00017236 0.0001669 0.00018098 | | | | | | | | | | 0.00020881 0.00017304 0.00018072 0.00017236 0.0001669 0.00018098 | | | | | 3.27135400e+03 3.27181620e+03 3.27605737e+03 3.27789380e+03 | | | | | | | | | | | | | 0.00019172 0.00017652 0.00016549 0.00021906 0.00022853 0.0001786 | | | | | | | | | | 0.00019172 0.00017652 0.00016549 0.00021906 0.00022853 0.0001786 | | | | | 3.27839380e+03 3.29614360e+03 3.33127657e+03 3.44105237e+03 | | | | | | | | | | | | | 0.00018281 0.00021391 0.00016893 0.00016621 0.00017098 0.00016548 | | | | | | | | | | 0.00018281 0.00021391 0.00016893 0.00016621 0.00017098 0.00016548 | | | | | 3.44461247e+03 3.50780473e+03 3.53047700e+03 3.70303903e+03 | | | | | | | | | | | | | 0.00018922 0.00017142 0.00017764 0.00016534 0.00017009 0.00017478 | | | | | | | | | | 0.00018922 0.00017142 0.00017764 0.00016534 0.00017009 0.00017478 | | | | | 3.70891697e+03 3.77949710e+03 3.85964867e+03 3.91479203e+03 | | | | | | | | | | | | | 0.00017044 0.00020593 0.00018119 0.00018153 0.00017829 0.00017439 | | | | | | | | | | 0.00017044 0.00020593 0.00018119 0.00018153 0.00017829 0.00017439 | | | | | 3.95843477e+03 3.96136523e+03 3.98141177e+03 3.98295023e+03 | | | | | | | | | | | | | 0.00023691 0.00017798 0.00016821 0.0002103 0.00016977 0.00018841 | | | | | | | | | | 0.00023691 0.00017798 0.00016821 0.0002103 0.00016977 0.00018841 | | | | | 4.00312307e+03 4.00379307e+03 4.00629407e+03 4.00779240e+03 | | | | | | | | | | | | | 0.00019102 0.0001787 0.00022766 0.00018956 0.00017718 0.00018092 | | | | | | | | | | 0.00019102 0.0001787 0.00022766 0.00018956 0.00017718 0.00018092 | | | | | 4.03416407e+03 4.04105260e+03 4.04658490e+03 4.04796910e+03 | | | | | | | | | | | | | 0.00016785 0.00017692 0.00017655 0.00017147 0.00019649 0.00018804 | | | | | | | | | | 0.00016785 0.00017692 0.00017655 0.00017147 0.00019649 0.00018804 | | | | | 4.06078403e+03 4.06175687e+03 4.06297263e+03 4.06595727e+03 | | | | | | | | | | | | | 0.000176 0.00019955 0.00017609 0.00018321 0.00016807 0.00019115 | | | | | | | | | | 0.000176 0.00019955 0.00017609 0.00018321 0.00016807 0.00019115 | | | | | 4.06753590e+03 4.10195333e+03 4.12225913e+03 4.13886977e+03 | | | | | | | | | | | | | 0.00016575 0.00016567 0.00022134 0.00027057 0.00017893 0.00016994 | | | | | | | | | | 0.00016575 0.00016567 0.00022134 0.00027057 0.00017893 0.00016994 | | | | | 4.14451150e+03 4.14570627e+03 4.15613140e+03 4.15698417e+03 | | | | | | | | | | | | | 0.00018647 0.00017711 0.00017535 0.00018941 0.00017954 0.0001845 | | | | | | | | | | 0.00018647 0.00017711 0.00017535 0.00018941 0.00017954 0.0001845 | | | | | 4.15745637e+03 4.15802600e+03 4.15865150e+03 4.16001690e+03 | | | | | | | | | | | | | 0.00019874 0.00018453 0.00017214 0.00017505 0.00018108 0.00016909 | | | | | | | | | | 0.00019874 0.00018453 0.00017214 0.00017505 0.00018108 0.00016909 | | | | | 4.16121103e+03 4.16623147e+03 4.16972737e+03 4.17580523e+03 | | | | | | | | | | | | | 0.00017405 0.00016654 0.00018489 0.00017906 0.00017702 0.00018694 | | | | | | | | | | 0.00017405 0.00016654 0.00018489 0.00017906 0.00017702 0.00018694 | | | | | 4.18053750e+03 4.19265073e+03 4.20319137e+03 4.20497980e+03 | | | | | | | | | | | | | 0.00016562 0.00018885 0.00017533 0.00018488] | | | | | | | | | | 0.00016562 0.00018885 0.00017533 0.00018488] | | | | | 4.20981620e+03 4.21801140e+03 4.24086563e+03 4.27133470e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.27829787e+03 4.27972480e+03 4.28175623e+03 4.31145593e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.31321843e+03 4.31380010e+03 4.33635120e+03 4.39187833e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.39586553e+03 4.42709677e+03 4.42791927e+03 4.42889893e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.42913487e+03 4.43059650e+03 4.43269933e+03 4.43734833e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.44594480e+03 4.45320670e+03 4.50156427e+03 4.51410553e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.52596233e+03 4.53713580e+03 4.54082017e+03 4.56210630e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.58856730e+03 4.61450923e+03 4.68068627e+03 4.68152283e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.68191117e+03 4.68202040e+03 4.68249737e+03 4.68358300e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.68727540e+03 4.68838623e+03 4.68986520e+03 4.69340210e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.69579457e+03 4.69837543e+03 4.70003177e+03 4.70216217e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.71068480e+03 4.71119830e+03 4.71259203e+03 4.71353037e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.71830693e+03 4.72500933e+03 4.72520693e+03 4.72785410e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.75477183e+03 4.78423573e+03 4.79557657e+03 4.79701650e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.81986893e+03 4.82318870e+03 4.82634597e+03 4.83272597e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.88208380e+03 4.88442477e+03 4.89177673e+03 4.89220653e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 4.92074810e+03 4.97544760e+03 5.00421063e+03 5.00426253e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.01475093e+03 5.01498230e+03 5.01538960e+03 5.01542843e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.01596743e+03 5.01648910e+03 5.01652907e+03 5.01732980e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.01874167e+03 5.02025217e+03 5.02297897e+03 5.02360487e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.02761127e+03 5.02778897e+03 5.02880223e+03 5.03054970e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.04421547e+03 5.05313337e+03 5.05746130e+03 5.05795433e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.05952560e+03 5.06374227e+03 5.06755533e+03 5.07873747e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.08227787e+03 5.08368817e+03 5.08807730e+03 5.10866150e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.11221523e+03 5.11666300e+03 5.11905597e+03 5.12879777e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.13023413e+03 5.13146840e+03 5.15065773e+03 5.16524393e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.16540767e+03 5.16584913e+03 5.17322130e+03 5.19487877e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.19491043e+03 5.19672020e+03 5.19776187e+03 5.20025047e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.21358653e+03 5.21374673e+03 5.21460957e+03 5.21490860e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.21538313e+03 5.21651280e+03 5.22656547e+03 5.22914390e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.23760140e+03 5.23905283e+03 5.25945437e+03 5.29977190e+03 | | | | | | | | | | | | | | | | | | | | | | | | | | | | 5.35360947e+03 5.35538160e+03 5.35823797e+03 5.35859823e+03] | | | | | | | | | | | | | | | | | | | | | | | | | | 2 | 664 | [6.79933333e-01 1.21883333e+00 1.59083333e+00 ... 5.35795553e+03 | 556 | 19.0886 | 4.36431 | -0.875027 | 1660 | 1 | MB | 0.000161646 | 164 | 1.56858 | probe00 | 1 | [0.00018863 0.00017188 0.000184 ... 0.00016821 0.00015934 0.00016874] | SCm | 0.000108756 | 19280.6 | 6.83513e-05 | 3245 | 0.605233 | 0.00799878 | 1.25356 | SCiw | [0.00018863 0.00017188 0.000184 ... 0.00016821 0.00015934 0.00016874] | 0.64432 | | | | 5.35799247e+03 5.35898827e+03] | | | | | | | | | | | | | | | | | | | | | | | | |
/opt/conda/lib/python3.10/site-packages/tabulate/__init__.py:108: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison or (len(row) >= 2 and row[1] == SEPARATING_LINE)
From it, we create a spikes_g object than can be rapidly indexed for analysis sake. It contains all the spike times and cluster number.
In [11]:
Copied!
spikes_g = dict()
input_list = good_units.spike_times.to_list()
spikes_g['times'] = np.concatenate(input_list).ravel()
spikes_list = list()
# We do not know how many spikes there is per row (this is not the spike count !)
# Loop over each unit
for index, spk in enumerate(input_list):
# Get the number of spikes per cluster
clu_no = good_units.unit_name[index]
clu = np.repeat(int(clu_no), spk.__len__())
# Append
spikes_list.append(clu)
spikes_g['clusters'] = np.concatenate(spikes_list)
# We now need to sort the spike times for analysis
index_sort = np.argsort(spikes_g['times'])
spikes_g['times'] = spikes_g['times'][index_sort]
spikes_g['clusters'] = spikes_g['clusters'][index_sort]
# Check the two attributes have the same size:
assert spikes_g["clusters"].__len__() == spikes_g["times"].__len__()
# Show example data for the first few spikes
print(spikes_g['clusters'][0:5])
print(spikes_g['times'][0:5])
spikes_g = dict()
input_list = good_units.spike_times.to_list()
spikes_g['times'] = np.concatenate(input_list).ravel()
spikes_list = list()
# We do not know how many spikes there is per row (this is not the spike count !)
# Loop over each unit
for index, spk in enumerate(input_list):
# Get the number of spikes per cluster
clu_no = good_units.unit_name[index]
clu = np.repeat(int(clu_no), spk.__len__())
# Append
spikes_list.append(clu)
spikes_g['clusters'] = np.concatenate(spikes_list)
# We now need to sort the spike times for analysis
index_sort = np.argsort(spikes_g['times'])
spikes_g['times'] = spikes_g['times'][index_sort]
spikes_g['clusters'] = spikes_g['clusters'][index_sort]
# Check the two attributes have the same size:
assert spikes_g["clusters"].__len__() == spikes_g["times"].__len__()
# Show example data for the first few spikes
print(spikes_g['clusters'][0:5])
print(spikes_g['times'][0:5])
[788 558 601 566 669] [0.02 0.02493333 0.02616667 0.03436667 0.03486667]
Channels¶
The channels dataframe contains the information for all the channels (e.g. AP and LF), for all the probes (here probe00 and probe01) acquired in a session.
In [12]:
Copied!
# ---------------------------------------------------
# Load channels data
channels_all = nwb.ec_electrodes.to_dataframe()
# Show names in dataframe
print(channels_all['channel_name'])
# ---------------------------------------------------
# Load channels data
channels_all = nwb.ec_electrodes.to_dataframe()
# Show names in dataframe
print(channels_all['channel_name'])
id
0 AP0
1 AP1
2 AP2
3 AP3
4 AP4
...
1531 LF379
1532 LF380
1533 LF381
1534 LF382
1535 LF383
Name: channel_name, Length: 1536, dtype: object
We will select only the AP channels for the probe of interest:
In [20]:
Copied!
# We select only the channels from the NeuropixelsShank00 -- Note this is different than probe_label
# And the channels in the AP band
channels = channels_all[(channels_all['group_name'] == 'NeuropixelsShank00') &
(channels_all['channel_name'].str.startswith('AP'))].copy()
# We add a column containing the atlas ID (Allen) to work with the IBL util plotting function
channels['atlas_id'] = ba.regions.acronym2id(channels['location']) # Convert from acronym to Allen ID
# We select only the channels from the NeuropixelsShank00 -- Note this is different than probe_label
# And the channels in the AP band
channels = channels_all[(channels_all['group_name'] == 'NeuropixelsShank00') &
(channels_all['channel_name'].str.startswith('AP'))].copy()
# We add a column containing the atlas ID (Allen) to work with the IBL util plotting function
channels['atlas_id'] = ba.regions.acronym2id(channels['location']) # Convert from acronym to Allen ID
Trials¶
The trials table is organised such that one row corresponds to one trial.
In [21]:
Copied!
# ---------------------------------------------------
# Load trial data
trials = nwb.trials.to_dataframe()
# Show first 3 rows of dataframe
print(trials.head(3).to_markdown())
# ---------------------------------------------------
# Load trial data
trials = nwb.trials.to_dataframe()
# Show first 3 rows of dataframe
print(trials.head(3).to_markdown())
| id | start_time | stop_time | choice | feedback_type | reward_volume | contrast_left | contrast_right | probability_left | feedback_time | response_time | stim_off_time | stim_on_time | go_cue_time | first_movement_time | |-----:|-------------:|------------:|---------:|----------------:|----------------:|----------------:|-----------------:|-------------------:|----------------:|----------------:|----------------:|---------------:|--------------:|----------------------:| | 0 | 15.6009 | 19.5275 | 1 | 1 | 1.5 | 1 | nan | 0.5 | 17.9712 | 17.9712 | 19.0273 | 16.9942 | 16.9952 | 17.7517 | | 1 | 19.9876 | 25.3734 | 1 | 1 | 1.5 | 0.125 | nan | 0.5 | 23.8064 | 23.8063 | 24.8734 | 23.4398 | 23.4408 | 23.6027 | | 2 | 25.8564 | 29.2062 | 1 | 1 | 1.5 | 1 | nan | 0.5 | 27.6402 | 27.6402 | 28.7061 | 27.3273 | 27.3281 | 27.4787 |
For the sake of this tutorial, we will remove any trials with NaNs in the first movement time:
In [22]:
Copied!
# Remove trials with nan values in first movement times
good_trials = trials[pd.notnull(trials.first_movement_time)]
# Remove trials with nan values in first movement times
good_trials = trials[pd.notnull(trials.first_movement_time)]
And we will convert each dataframe column into numpy arrays:
In [23]:
Copied!
# ---------------------------------------------------
# Convert to numpy for analysis
events = good_trials.first_movement_time.to_numpy()
contrast_R = good_trials.contrast_right.to_numpy()
contrast_L = good_trials.contrast_left.to_numpy()
choice = good_trials.choice.to_numpy()
block = good_trials.probability_left.to_numpy()
# N trial count
num_trial = len(events)
# Find "trials" that go in one direction and the other direction
# Note: This is not a pure indexing on the *task trials* as we removed trials with nan values previously
indx_choice_a = np.where(choice == -1)[0]
indx_choice_b = np.where(choice == 1)[0]
# ---------------------------------------------------
# Convert to numpy for analysis
events = good_trials.first_movement_time.to_numpy()
contrast_R = good_trials.contrast_right.to_numpy()
contrast_L = good_trials.contrast_left.to_numpy()
choice = good_trials.choice.to_numpy()
block = good_trials.probability_left.to_numpy()
# N trial count
num_trial = len(events)
# Find "trials" that go in one direction and the other direction
# Note: This is not a pure indexing on the *task trials* as we removed trials with nan values previously
indx_choice_a = np.where(choice == -1)[0]
indx_choice_b = np.where(choice == 1)[0]
Wheel¶
The data saved in NWB is the raw wheel data. In order to get the wheel position and velocity, we need to filter this raw data with the ibl util functions.
In [24]:
Copied!
# ---------------------------------------------------
# Load the wheel data
wheel_raw_position = nwb.processing["behavior"].containers["CompassDirection"].spatial_series["WheelPositionSeries"].data
wheel_raw_timestamps = nwb.processing["behavior"].containers["CompassDirection"].spatial_series["WheelPositionSeries"].timestamps
wheel = pd.DataFrame(columns=['timestamps', 'position', 'velocity', 'acceleration'])
wheel['position'], wheel['timestamps'] = interpolate_position(wheel_raw_timestamps, wheel_raw_position, freq=1000)
wheel['velocity'], wheel['acceleration'] = velocity_filtered(wheel['position'], fs=1000, corner_frequency=20, order=8)
wheel = wheel.apply(np.float32)
# Show first 3 rows of dataframe
print(wheel.head(3).to_markdown())
# ---------------------------------------------------
# Load the wheel data
wheel_raw_position = nwb.processing["behavior"].containers["CompassDirection"].spatial_series["WheelPositionSeries"].data
wheel_raw_timestamps = nwb.processing["behavior"].containers["CompassDirection"].spatial_series["WheelPositionSeries"].timestamps
wheel = pd.DataFrame(columns=['timestamps', 'position', 'velocity', 'acceleration'])
wheel['position'], wheel['timestamps'] = interpolate_position(wheel_raw_timestamps, wheel_raw_position, freq=1000)
wheel['velocity'], wheel['acceleration'] = velocity_filtered(wheel['position'], fs=1000, corner_frequency=20, order=8)
wheel = wheel.apply(np.float32)
# Show first 3 rows of dataframe
print(wheel.head(3).to_markdown())
| | timestamps | position | velocity | acceleration | |---:|-------------:|-----------:|-----------:|---------------:| | 0 | 11.0097 | 0.00153398 | 0 | 0 | | 1 | 11.0107 | 0.0014506 | -0.0918616 | -91.8616 | | 2 | 11.0117 | 0.00136723 | -0.0917454 | 0.11619 |
Note: You can access the wheel data through the nwb.modules['behavior'] module, for example
velocity_df = nwb.modules['behavior'].data_interfaces['WheelVelocity'].to_dataframe()
🔴 however this data has been processed differently than what is done at the IBL, and therefore we recommend you to apply the filtering functions above on the raw wheel data to obtain wheel position, speed, velocity and acceleration.
Encoding¶
Introduction¶
Here we want to assess whether single units encode the choice of the animal (i.e. turning the wheel left or right).
To not be confounded by the encoding of the movement itself, we will use the neural activity solely in a narrow time-window prior to the first movement onset as detected on the wheel trace.
We will compute in this time window the modulation index, i.e. the change in firing rate between the trials where the animal turned the wheel left or right.
We will assess the significance of such modulation by doing complex permutation tests.
Analysis¶
Here we use a time window of 100ms prior to the first movement onset:
In [25]:
Copied!
# ---------------------------------------------------
# Select a time window of interest
time_window = np.array([-0.1, 0.0]) # 100 ms before the event
# ---------------------------------------------------
# Select a time window of interest
time_window = np.array([-0.1, 0.0]) # 100 ms before the event
Let's display the wheel position and speed for the first 4 trials, to see where the time window is compared to it:
In [26]:
Copied!
# ---------------------------------------------------
# Plot the whole wheel trace for the first N events in the task
n_ev_plt = 4
indx_plot = np.where((wheel.timestamps > events[0] - 1.0) & (wheel.timestamps < events[n_ev_plt-1] + 1.0))[0]
# Positive wheel speed indicate a CCW (left stim) move
# -1 (turn CCW), +1 (turn CW), or 0 (nogo)
choice_text = np.repeat('CW (left stim -> center)', n_ev_plt)
choice_text[np.where(choice[0:n_ev_plt] == -1)] = 'CCW (right stim -> center)'
# We represent left choice in red / right choice in blue
choice_color = np.repeat('#0000FF', n_ev_plt)
choice_color[np.where(choice[0:n_ev_plt] == -1)] = '#FF0000'
fig, ax = plt.subplots(1)
fig.set_size_inches(18.5, 7.5)
plt.plot(wheel.timestamps[indx_plot], wheel.velocity[indx_plot], 'k', label='wheel speed')
plt.plot(wheel.timestamps[indx_plot], wheel.position[indx_plot], label='wheel position')
plt.scatter(events[0:n_ev_plt], 4*np.ones(events[0:n_ev_plt].shape), color=choice_color, marker='v')
for i_ev_plt in range(0, n_ev_plt):
if not np.isnan(contrast_R[i_ev_plt]):
contrast_text = f'Contrast right: {contrast_R[i_ev_plt]*100}%'
else:
contrast_text = f'Contrast left: {contrast_L[i_ev_plt] * 100}%'
txt = f'Choice {choice_text[i_ev_plt]} \n{contrast_text}'
ax.annotate(txt, (events[i_ev_plt]-1.5, 4.5))
# Highlight the part where the ephys data is taken with rectangles
for event in events[0:n_ev_plt]:
ax.add_patch(Rectangle((event+time_window[0], -4), width=time_window[1]-time_window[0],
height=8, alpha=0.5, color='#00FFAA', edgecolor=None))
ax.legend(loc=4)
# ---------------------------------------------------
# Plot the whole wheel trace for the first N events in the task
n_ev_plt = 4
indx_plot = np.where((wheel.timestamps > events[0] - 1.0) & (wheel.timestamps < events[n_ev_plt-1] + 1.0))[0]
# Positive wheel speed indicate a CCW (left stim) move
# -1 (turn CCW), +1 (turn CW), or 0 (nogo)
choice_text = np.repeat('CW (left stim -> center)', n_ev_plt)
choice_text[np.where(choice[0:n_ev_plt] == -1)] = 'CCW (right stim -> center)'
# We represent left choice in red / right choice in blue
choice_color = np.repeat('#0000FF', n_ev_plt)
choice_color[np.where(choice[0:n_ev_plt] == -1)] = '#FF0000'
fig, ax = plt.subplots(1)
fig.set_size_inches(18.5, 7.5)
plt.plot(wheel.timestamps[indx_plot], wheel.velocity[indx_plot], 'k', label='wheel speed')
plt.plot(wheel.timestamps[indx_plot], wheel.position[indx_plot], label='wheel position')
plt.scatter(events[0:n_ev_plt], 4*np.ones(events[0:n_ev_plt].shape), color=choice_color, marker='v')
for i_ev_plt in range(0, n_ev_plt):
if not np.isnan(contrast_R[i_ev_plt]):
contrast_text = f'Contrast right: {contrast_R[i_ev_plt]*100}%'
else:
contrast_text = f'Contrast left: {contrast_L[i_ev_plt] * 100}%'
txt = f'Choice {choice_text[i_ev_plt]} \n{contrast_text}'
ax.annotate(txt, (events[i_ev_plt]-1.5, 4.5))
# Highlight the part where the ephys data is taken with rectangles
for event in events[0:n_ev_plt]:
ax.add_patch(Rectangle((event+time_window[0], -4), width=time_window[1]-time_window[0],
height=8, alpha=0.5, color='#00FFAA', edgecolor=None))
ax.legend(loc=4)
Out[26]:
<matplotlib.legend.Legend at 0x7fe5dbd9f220>
Now we compute the spike rate around our event of interest (first movement), and then the modulation index indicating the magnitude of the change in firing rate between the two conditions (choice left or right). A high absolute value of the modulation index indicates a strong change between the firing in the two conditions, and therefore a potential encoding.
We also perform a significance test to obtain a p-value showing whether the activity is different for left versus right choices :
In [27]:
Copied!
# ---------------------------------------------------
# Compute spike rate around event
events_tw = np.array([events+time_window[0], events+time_window[1]]).T
# Compute count (for all clusters of interest) (THIS CAN TAKE A WHILE)
spike_count, cluster_id = get_spike_counts_in_bins(spikes_g['times'], spikes_g['clusters'], events_tw)
# Note: the cluster_id returned are stricty increasing, as we have in the dataframe:
print(f'the cluster_id returned :\n\n {cluster_id} \n\n is in the same order as in the dataframe :\n')
print(good_units.unit_name.to_numpy().astype('int'))
# Validate
np.testing.assert_equal(cluster_id, good_units.unit_name.to_numpy().astype('int'))
# Compute rate (for all clusters of interest)
spike_rate = np.zeros((num_neuron, num_trial))
spike_rate = spike_count / (time_window[1] - time_window[0])
# ---------------------------------------------------
# Compute the modulation index
modulation_index = np.divide((np.mean(spike_rate[:, indx_choice_a], 1) - np.mean(spike_rate[:, indx_choice_b], 1)),
(np.mean(spike_rate[:, indx_choice_a], 1) + np.mean(spike_rate[:, indx_choice_b], 1)))
# ---------------------------------------------------
# Compute spike rate around event
events_tw = np.array([events+time_window[0], events+time_window[1]]).T
# Compute count (for all clusters of interest) (THIS CAN TAKE A WHILE)
spike_count, cluster_id = get_spike_counts_in_bins(spikes_g['times'], spikes_g['clusters'], events_tw)
# Note: the cluster_id returned are stricty increasing, as we have in the dataframe:
print(f'the cluster_id returned :\n\n {cluster_id} \n\n is in the same order as in the dataframe :\n')
print(good_units.unit_name.to_numpy().astype('int'))
# Validate
np.testing.assert_equal(cluster_id, good_units.unit_name.to_numpy().astype('int'))
# Compute rate (for all clusters of interest)
spike_rate = np.zeros((num_neuron, num_trial))
spike_rate = spike_count / (time_window[1] - time_window[0])
# ---------------------------------------------------
# Compute the modulation index
modulation_index = np.divide((np.mean(spike_rate[:, indx_choice_a], 1) - np.mean(spike_rate[:, indx_choice_b], 1)),
(np.mean(spike_rate[:, indx_choice_a], 1) + np.mean(spike_rate[:, indx_choice_b], 1)))
the cluster_id returned : [545 548 556 558 565 566 581 582 584 587 601 603 604 610 613 640 659 669 682 683 684 686 697 706 709 721 724 740 746 754 764 770 775 782 785 788 799 800 807 808 810 811 812 815 820 832 834 836 846 849 858] is in the same order as in the dataframe : [545 548 556 558 565 566 581 582 584 587 601 603 604 610 613 640 659 669 682 683 684 686 697 706 709 721 724 740 746 754 764 770 775 782 785 788 799 800 807 808 810 811 812 815 820 832 834 836 846 849 858]
In [28]:
Copied!
# ---------------------------------------------------
# Compute p-value to assess significance of encoding (THIS CAN TAKE A WHILE)
p_1 = get_choice_time_shuffle(spike_rate, contrast_L, contrast_R, block, choice, 1000)
# Find those that are significant to 0.01
p_sig = np.zeros(p_1.shape)
p_sig[np.where(p_1 < 0.01)] = 1
# ---------------------------------------------------
# Compute p-value to assess significance of encoding (THIS CAN TAKE A WHILE)
p_1 = get_choice_time_shuffle(spike_rate, contrast_L, contrast_R, block, choice, 1000)
# Find those that are significant to 0.01
p_sig = np.zeros(p_1.shape)
p_sig[np.where(p_1 < 0.01)] = 1
Now we plot the modulation index value for each neuron along the probe depth, colored according to the level of significance.
We can see that a few units are significantly modulated.
In [29]:
Copied!
# ---------------------------------------------------
# Plot the modulation index for each unit across brain regions
fig, axs = plt.subplots(1, 2)
fig.set_size_inches(6.5, 5.5)
# Plot the brain regions through which the channels on the insertion pass
plot_brain_regions(channels.atlas_id.to_numpy(), channel_depths=channels.rel_y.to_numpy(), ax=axs[0])
# Plot scatter plot of cluster depths vs cluster modulation index
axs[1].scatter(modulation_index, good_units.mean_relative_depth, c=p_sig) # color by significance
axs[1].set_xlabel('M.I.')
axs[1].get_yaxis().set_visible(False)
# ---------------------------------------------------
# Plot the modulation index for each unit across brain regions
fig, axs = plt.subplots(1, 2)
fig.set_size_inches(6.5, 5.5)
# Plot the brain regions through which the channels on the insertion pass
plot_brain_regions(channels.atlas_id.to_numpy(), channel_depths=channels.rel_y.to_numpy(), ax=axs[0])
# Plot scatter plot of cluster depths vs cluster modulation index
axs[1].scatter(modulation_index, good_units.mean_relative_depth, c=p_sig) # color by significance
axs[1].set_xlabel('M.I.')
axs[1].get_yaxis().set_visible(False)
Note to go further:
- This analysis can be done using any time window, or any events. As such, you can try to see if units are also selective to left / right conditions after the first move, or after the stimulus onset for example.
Manifold¶
Introduction¶
In the above example (Encoding), we looked at whether units encoded significantly the choice (left or right) using their firing rate in a time window before the first movement.
Now, we will study when the neural activity representing left or right choices differs before the first movement.
We will use a dimensionality approach (Manifold) to measure the distance between the left and right choice representation during the time window of interest.
The manifold analysis steps are: (see figure below)
- a) For each unit, take the spikes before the first movement in a time window, and create a PSTH for left and right choices
- b) Stack these average PSTHs across units so as to form a high-dimentional representation of the choice (N dimension = N units)
- c) Reduce these PSTHs to a lower-dimention state space, e.g. through PCA, and plot the trajectories in that space (using e.g. the first 3 PCs)
- d) Compute the distance between the reduced trajectories at each timepoint, and assess the significance of the distance (which would mean, whether the units represent the choice left/right above chance)
- e) Find the time point for which the distance is at a percentage (e.g. 70%) of the maximal distance (
latency)- Note: do not use the absolute max otherwise most latencies will be at the first movement time
In (c-e) we show you the results for all units in the region GRN (gigantocellular reticular nucleus). To do this, we had to load units from multiple probe insertions targeting the GRN.
In (f-g) we show you the result for all regions of the Brainwide map. The code to reproduce this figure written by Michael Schartner is here.
Note: This analysis can be done using single-units for a given region from a single PID as is done in the example below, however it works best if you stack as many PSTHs as possible. As such, it is better to first load the spikes data from many recordings spanning a region of interest at once.
Analysis¶
Here we will do a much simplified version of the analysis, without any control to show you how to start simply.
We will use PCA, and focus only on units in the motor-area of the Superior Colliculus (SCm).
Select the units in a given brain region (here SCm):
In [30]:
Copied!
# ---------------------------------------------------
# Select units in a given brain region
# Print each unit's brain region label, for the first 10 units as example
print(good_units['beryl_location'].head(10).to_markdown())
# Use the Beryl brain region parcellation to select those in SCm
units_SCm = good_units[good_units['beryl_location'] == 'SCm']
# Take the index of units in SCm
cluster_SCm_IDs = units_SCm.unit_name.to_numpy().astype('int')
nunit = len(cluster_SCm_IDs)
print(f'\n We will be working with {nunit} SCm units')
# ---------------------------------------------------
# Select units in a given brain region
# Print each unit's brain region label, for the first 10 units as example
print(good_units['beryl_location'].head(10).to_markdown())
# Use the Beryl brain region parcellation to select those in SCm
units_SCm = good_units[good_units['beryl_location'] == 'SCm']
# Take the index of units in SCm
cluster_SCm_IDs = units_SCm.unit_name.to_numpy().astype('int')
nunit = len(cluster_SCm_IDs)
print(f'\n We will be working with {nunit} SCm units')
| | beryl_location | |---:|:-----------------| | 0 | SCm | | 1 | SCm | | 2 | SCm | | 3 | RN | | 4 | RSPd | | 5 | MRN | | 6 | MRN | | 7 | MRN | | 8 | SCm | | 9 | SCm | We will be working with 21 SCm units
(a-b) Create stacked PSTHs (using all cells):
In [31]:
Copied!
# ---------------------------------------------------
# Create PSTHs
binsize = 0.01 # bin size [sec] for neural binning
time_window = [-0.150, 0]
for count, clu_id in enumerate(cluster_SCm_IDs):
# Find spikes for this cluster
spk_indx = np.where(np.isin(spikes_g['clusters'], clu_id))
spikes_unit = {key: val[spk_indx] for key, val in spikes_g.items()}
# Compute raster
raster, timestamps = get_event_aligned_raster(spikes_unit['times'], events, tbin=binsize, values=None,
epoch=time_window, bin=True)
# Compute PSTH (return only the mean)
psth_a, _ = get_psth(raster, trial_ids=indx_choice_a)
psth_b, _ = get_psth(raster, trial_ids=indx_choice_b)
# ------- Stack PSTHs -------
# Init ; Here we create a M = n condition x n unit x n time bin and will concatenate it later
if count == 0:
nbin = len(timestamps)
stack_psth = np.empty((2, nunit, nbin))
stack_psth[0, count, :] = psth_a
stack_psth[1, count, :] = psth_b
# ---------------------------------------------------
# Create PSTHs
binsize = 0.01 # bin size [sec] for neural binning
time_window = [-0.150, 0]
for count, clu_id in enumerate(cluster_SCm_IDs):
# Find spikes for this cluster
spk_indx = np.where(np.isin(spikes_g['clusters'], clu_id))
spikes_unit = {key: val[spk_indx] for key, val in spikes_g.items()}
# Compute raster
raster, timestamps = get_event_aligned_raster(spikes_unit['times'], events, tbin=binsize, values=None,
epoch=time_window, bin=True)
# Compute PSTH (return only the mean)
psth_a, _ = get_psth(raster, trial_ids=indx_choice_a)
psth_b, _ = get_psth(raster, trial_ids=indx_choice_b)
# ------- Stack PSTHs -------
# Init ; Here we create a M = n condition x n unit x n time bin and will concatenate it later
if count == 0:
nbin = len(timestamps)
stack_psth = np.empty((2, nunit, nbin))
stack_psth[0, count, :] = psth_a
stack_psth[1, count, :] = psth_b
Plot the stacked PSTHs:
In [32]:
Copied!
# ---------------------------------------------------
# Plot stacked PSTHs
fig, ax = plt.subplots(2)
fig.set_size_inches(6.5, 5.5)
ax[0].imshow(stack_psth[0, :, :], vmax=stack_psth.max(), vmin=stack_psth.min())
ax[1].imshow(stack_psth[1, :, :], vmax=stack_psth.max(), vmin=stack_psth.min())
ax[0].set_xlabel('')
ax[1].set_xlabel('time bin')
ax[0].set_ylabel('PSTH CCW (right) choice')
ax[1].set_ylabel('PSTH CW (left) choice')
ytick_loc = [0]
ax[0].set_yticks(ytick_loc)
ax[1].set_yticks(ytick_loc)
ax[0].set_yticklabels(['unit #1'])
ax[1].set_yticklabels(['unit #1'])
# ---------------------------------------------------
# Plot stacked PSTHs
fig, ax = plt.subplots(2)
fig.set_size_inches(6.5, 5.5)
ax[0].imshow(stack_psth[0, :, :], vmax=stack_psth.max(), vmin=stack_psth.min())
ax[1].imshow(stack_psth[1, :, :], vmax=stack_psth.max(), vmin=stack_psth.min())
ax[0].set_xlabel('')
ax[1].set_xlabel('time bin')
ax[0].set_ylabel('PSTH CCW (right) choice')
ax[1].set_ylabel('PSTH CW (left) choice')
ytick_loc = [0]
ax[0].set_yticks(ytick_loc)
ax[1].set_yticks(ytick_loc)
ax[0].set_yticklabels(['unit #1'])
ax[1].set_yticklabels(['unit #1'])
Out[32]:
[Text(0, 0, 'unit #1')]
(c) Generate state space:
In [33]:
Copied!
# ---------------------------------------------------
# Perform PCA
pca = PCA(n_components=2)
# Use concatenate to place the PSTHs of the two conditions one after the other in time
# Transpose before doing the PCA, so it's done on the unit axis
trajs = pca.fit_transform(np.concatenate(stack_psth, axis=1).T).T
# ---------------------------------------------------
# Take trajectory per condition
traj_a = trajs[:, 0:nbin]
traj_b = trajs[:, nbin:]
# ---------------------------------------------------
# Perform PCA
pca = PCA(n_components=2)
# Use concatenate to place the PSTHs of the two conditions one after the other in time
# Transpose before doing the PCA, so it's done on the unit axis
trajs = pca.fit_transform(np.concatenate(stack_psth, axis=1).T).T
# ---------------------------------------------------
# Take trajectory per condition
traj_a = trajs[:, 0:nbin]
traj_b = trajs[:, nbin:]
In [34]:
Copied!
# ---------------------------------------------------
# Plot the trajectories
fig, ax = plt.subplots(1)
fig.set_size_inches(6.5, 5.5)
ax.plot(traj_a[0, :], traj_a[1, :], color='b', label='CCW (right)', marker='.')
ax.plot(traj_b[0, :], traj_b[1, :], color='r', label='CW (left)', marker='.')
# Mark the start
ax.plot(traj_a[0, 0], traj_a[1, 0], color='b', marker='*')
ax.plot(traj_b[0, 0], traj_b[1, 0], color='r', marker='*')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.legend(loc=2)
# ---------------------------------------------------
# Plot the trajectories
fig, ax = plt.subplots(1)
fig.set_size_inches(6.5, 5.5)
ax.plot(traj_a[0, :], traj_a[1, :], color='b', label='CCW (right)', marker='.')
ax.plot(traj_b[0, :], traj_b[1, :], color='r', label='CW (left)', marker='.')
# Mark the start
ax.plot(traj_a[0, 0], traj_a[1, 0], color='b', marker='*')
ax.plot(traj_b[0, 0], traj_b[1, 0], color='r', marker='*')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.legend(loc=2)
Out[34]:
<matplotlib.legend.Legend at 0x7fe5db33ac50>
Let's now compute and plot the Euclidean distance between the two trajectories:
In [35]:
Copied!
# ---------------------------------------------------
# Compute the Euclidean distance
euc_dist = np.sqrt(np.sum((traj_a - traj_b) ** 2, axis=0))
# ---------------------------------------------------
# Plot the euclidean distance
fig, ax = plt.subplots(1)
fig.set_size_inches(6.5, 5.5)
ax.plot(timestamps, euc_dist, color='k', label='SCm')
ax.set_xlabel('time')
ax.set_ylabel('Euclidian distance')
# ---------------------------------------------------
# Compute the Euclidean distance
euc_dist = np.sqrt(np.sum((traj_a - traj_b) ** 2, axis=0))
# ---------------------------------------------------
# Plot the euclidean distance
fig, ax = plt.subplots(1)
fig.set_size_inches(6.5, 5.5)
ax.plot(timestamps, euc_dist, color='k', label='SCm')
ax.set_xlabel('time')
ax.set_ylabel('Euclidian distance')
Out[35]:
Text(0, 0.5, 'Euclidian distance')
Note that this is rather wobbly -- using more cells from other PIDs would help smoothen the curves.
Decoding¶
In this example, we will want to decode the choice for a given trial, using the activity of several units (you can select units from a chosen brain region, here for simplicity we use all units).
As the choice is a binary variable (+1, -1) we will use logistic regression in a high dimension space (defined by the number of units).
We will devise training and testing sets as is classically done to perform accuracy assessment in machine learning.
Useful reference:
- To prepare test and training sets, use the built-in
scikit-learn train_test_splitfunction.
In [36]:
Copied!
spike_count.shape
index_trials = np.where((choice == 1) | (choice == -1))[0]
len(index_trials)
spike_count[:, index_trials].shape
spike_count.shape
index_trials = np.where((choice == 1) | (choice == -1))[0]
len(index_trials)
spike_count[:, index_trials].shape
Out[36]:
(51, 691)
In [37]:
Copied!
# ---------------------------------------------------
# Check that there are indeed only 2 values in choice
if len(np.unique(choice)) != 2:
# Find trials that have a choice value of +1 or -1
index_trials = np.where((choice == 1) | (choice == -1))[0]
# Reassign to only those trials
choice = choice[index_trials]
spike_count = spike_count[:, index_trials]
# Ensure there are indeed two values now
assert len(np.unique(choice)) == 2
# ---------------------------------------------------
# Use the spike count in the time window prior to the movement onset as predictors.
# We use the same time window as defined earlier (Encoding)
X = spike_count.T # shape of spike_count : n units x n trials in set -> transpose to fit model
y = choice
# ---------------------------------------------------
# Split trials into test and training sets for the logistic regression.
# Take first half of trials for training, second half for testing (test_size = 0.5)
# Fix random seed to repeat results across runs
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)
# ---------------------------------------------------
# See how many trials there is for choice +1 ; -1
print(f'Training set : {len(np.where(y_train == 1)[0])} choice +1 / {len(y_train)} trials')
print(f'Testing set : {len(np.where(y_test == 1)[0])} choice +1 / {len(y_test)} trials')
# ---------------------------------------------------
# Check that there are indeed only 2 values in choice
if len(np.unique(choice)) != 2:
# Find trials that have a choice value of +1 or -1
index_trials = np.where((choice == 1) | (choice == -1))[0]
# Reassign to only those trials
choice = choice[index_trials]
spike_count = spike_count[:, index_trials]
# Ensure there are indeed two values now
assert len(np.unique(choice)) == 2
# ---------------------------------------------------
# Use the spike count in the time window prior to the movement onset as predictors.
# We use the same time window as defined earlier (Encoding)
X = spike_count.T # shape of spike_count : n units x n trials in set -> transpose to fit model
y = choice
# ---------------------------------------------------
# Split trials into test and training sets for the logistic regression.
# Take first half of trials for training, second half for testing (test_size = 0.5)
# Fix random seed to repeat results across runs
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)
# ---------------------------------------------------
# See how many trials there is for choice +1 ; -1
print(f'Training set : {len(np.where(y_train == 1)[0])} choice +1 / {len(y_train)} trials')
print(f'Testing set : {len(np.where(y_test == 1)[0])} choice +1 / {len(y_test)} trials')
Training set : 272 choice +1 / 345 trials Testing set : 266 choice +1 / 346 trials
In [38]:
Copied!
# ---------------------------------------------------
# Fit the logistic regression model
clf = LogisticRegression(random_state=0).fit(X_train, y_train)
# Use the test set to assess the model accuracy
choice_predicted = clf.predict(X_test)
n_trial_correct = len(np.where(y_test == choice_predicted)[0])
print(f'Accuracy : {n_trial_correct} trials correctly predicted / {len(y_test)} trials')
# ---------------------------------------------------
# Fit the logistic regression model
clf = LogisticRegression(random_state=0).fit(X_train, y_train)
# Use the test set to assess the model accuracy
choice_predicted = clf.predict(X_test)
n_trial_correct = len(np.where(y_test == choice_predicted)[0])
print(f'Accuracy : {n_trial_correct} trials correctly predicted / {len(y_test)} trials')
Accuracy : 286 trials correctly predicted / 346 trials
/opt/conda/lib/python3.10/site-packages/sklearn/linear_model/_logistic.py:458: ConvergenceWarning: lbfgs failed to converge (status=1):
STOP: TOTAL NO. of ITERATIONS REACHED LIMIT.
Increase the number of iterations (max_iter) or scale the data as shown in:
https://scikit-learn.org/stable/modules/preprocessing.html
Please also refer to the documentation for alternative solver options:
https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression
n_iter_i = _check_optimize_result(
Note to go further:
Here we split the trials half/half for training and testing set. To get a better sense for accuracy, you could train using all but 1 trial (leave-one-out), and run the same analysis as above in a loop over all trials - this way, you would use as much information as possible to fit the logistic regression.
References
- A much more thorough analysis example can be found here; here is the link to its most inner fitting function; code written by Matt Whiteway and Brandon Benson.
More ideas¶
📗 We studied choice - what about other binary variables?¶
Are you interested to know how units encode other variables, such as the reward or the visual stimulus ?
You can run similar analysis as the Encoding, Manifold and Decoding, but for other types of binary variables present in the task (printed below), such as :
- the
visual stimulus side(presented on the left or right side of the screen) - the
feedback type(reward or noise burst) - the
movement side(CCW or CW)
Respectively, you can use these variables to create a time window of interest:
- (0, +0.1) second from
stim_on_time - (0, +0.1) second from
feedback_time - (0, +0.1) second from
first_movement_time(note that this time the window is after the first movement)
In [ ]:
Copied!
trials.keys()
trials.keys()
We show you below interesting brain regions to look at using the example of the manifold analysis:
- In (a) is shown the trajectory of the brain region
GRNfor the choice, which was our example variable in this tutorial. - In (b) are shown two examples regions,
VISpandIRNthat are encoding thestimulus sideandfeedback typerespectively.- Note that the time window taken here is positive (i.e. after the stimulus or feedback onset), and much longer for the case of the feedback.
- Note that IRN presents an oscillatory pattern which may be linked to licking.
Useful links:
- Documentation on how to find a PID going through a brain region of interest ; you can also use this query, which returns a list of PIDs:
# TODO WITH NWB SEARCH
The behavior module of NWB contains the lick time:
In [ ]:
Copied!
print(nwb.modules['behavior'])
licks = nwb.modules['behavior'].data_interfaces['LickTimes'].to_dataframe()
print(nwb.modules['behavior'])
licks = nwb.modules['behavior'].data_interfaces['LickTimes'].to_dataframe()
📙 Delineating mixed-responses using General Linear Model (GLM)¶
Introduction¶
In the Encoding example above, we showed you how to probe whether units were significantly modulated between two conditions of a given event in a specific time window (e.g. event: choice, condition: left vs right, time window: 100ms prior to the first movement).
However, units can be modulated by several events in the task (e.g. visual stimuli and choice). In our task these events often occur very close temporally, as the reaction time of the mouse can be as quick as a few hundreds of milliseconds.
To disambiguate the contribution of each task variables, you can fit response kernels mimicking the time course neural responses using a General Linera Model (GLM) to the real neural responses.
The weight of each kernel indicates how much a unit is responsive to a given task event.
Example figure and code¶
In the plot below, we show an example fit for several kernels to one unit. You can see this unit is mostly encoding the right stimulus onset, as most other kernels are flat.
To make these plots, we use code from the neuroencoding repository, created by Berk Gercek.
Useful links from within this repository:
- Example code to run a GLM fit for a single session
- Code for the design matrix creation (this is the hard part)
- Use the Class GLMPredictor to create plots
How to install the repository
🐍 Note: This may not work in Colab, install this on your machine instead 🪛
In [ ]:
Copied!
!git clone https://github.com/berkgercek/neurencoding.git
!git clone https://github.com/berkgercek/neurencoding.git
In [ ]:
Copied!
%cd /content/neurencoding
!pip install -e .
%cd /content/neurencoding
!pip install -e .
📘 Correlate neural activity with continuous variables instead of binary¶
Some variables in the task are binary (e.g. choice, stimulus side), but some others are continuous.
We list below a few continuous variables:
- the
wheel speedandposition - the
paw position(loading documentation here) TODO CHANGE FOR NWB DOC - the
motion energy(loading documentation here)
Using e.g. linear regression, you can try to decode the values of these continuous variables using neural data.
A question that such analysis would answer would be:
- Do neurons encode the vigor of movement?
📕 Model the behavior and see if there are neural correlates¶
How would you answer the questions:
- Do animal adopt a certain body posture, that is predictive of their decision?
- Where are these postures encoded in the brain ?
Use the output from the video segmentation algorithm (loading documentation here), e.g. the paw position (paw_R) to see if there are any stereotypical movements.
Use an unsupervised classifier method to segregate different kinds of movement the animals makes, and assess whether there are particular patterns of occurrence.
Assess whether these stereotypical movements can be decoded from the neural activity using the techniques of this tutorial (Encoding or Decoding).
Read further on how to access and load the data¶
- Cosyne tutorial on data access
- Loading data examples TODO change this for NWB page